Simulations of complementarity systems using time-stepping
One of the useful outcomes of the theory of complementarity systems is a new family of methods for numerical simulation of discontinuous systems. Here we will demonstrate the essence by introducing the method of time-stepping. And we do it by means of an example.
Example 1 (Simulation using time-stepping) Consider the following discontinuous dynamical system in \mathbb R^2:
\begin{aligned}
\dot x_1 &= -\operatorname{sign} x_1 + 2 \operatorname{sign} x_2\\
\dot x_2 &= -2\operatorname{sign} x_1 -\operatorname{sign} x_2.
\end{aligned}
Precompiling CairoMakie...
612.6 ms ✓ OffsetArrays
507.1 ms ✓ OffsetArrays → OffsetArraysAdaptExt
459.7 ms ✓ StackViews
529.4 ms ✓ PaddedViews
592.3 ms ✓ MosaicViews
1433.7 ms ✓ Interpolations
907.9 ms ✓ Interpolations → InterpolationsUnitfulExt
1079.0 ms ✓ KernelDensity
8627.1 ms ✓ ImageCore
1617.2 ms ✓ ImageBase
2200.1 ms ✓ Sixel
2442.6 ms ✓ JpegTurbo
2783.7 ms ✓ PNGFiles
1482.5 ms ✓ WebP
1199.4 ms ✓ ImageAxes
778.6 ms ✓ ImageMetadata
1152.1 ms ✓ Netpbm
27596.9 ms ✓ TiffImages
65816.2 ms ✓ Makie
23175.5 ms ✓ CairoMakie
20 dependencies successfully precompiled in 118 seconds. 255 already precompiled.
Figure 1: State portrait of the discontinuous system
One particular (vector) state trajectory is in Fig. 2.
Precompiling DifferentialEquations...
320.8 ms ✓ StaticArrayInterface → StaticArrayInterfaceOffsetArraysExt
636.8 ms ✓ Sparspak
12044.4 ms ✓ LoopVectorization
781.0 ms ✓ LoopVectorization → SpecialFunctionsExt
855.3 ms ✓ LoopVectorization → ForwardDiffExt
853.0 ms ✓ TriangularSolve
4506.7 ms ✓ RecursiveFactorization
13509.2 ms ✓ LinearSolve
2045.7 ms ✓ LinearSolve → LinearSolveRecursiveArrayToolsExt
2075.7 ms ✓ LinearSolve → LinearSolveBandedMatricesExt
2089.7 ms ✓ LinearSolve → LinearSolveEnzymeExt
2161.6 ms ✓ NonlinearSolveBase → NonlinearSolveBaseLinearSolveExt
2169.1 ms ✓ LinearSolve → LinearSolveFastAlmostBandedMatricesExt
2523.2 ms ✓ OrdinaryDiffEqDifferentiation
3907.6 ms ✓ OrdinaryDiffEqExtrapolation
11369.7 ms ✓ NonlinearSolveFirstOrder
11997.4 ms ✓ OrdinaryDiffEqRosenbrock
3765.9 ms ✓ BoundaryValueDiffEqCore
25735.2 ms ✓ NonlinearSolve
3877.0 ms ✓ NonlinearSolve → NonlinearSolveBandedMatricesExt
3878.6 ms ✓ NonlinearSolve → NonlinearSolveNLsolveExt
4730.2 ms ✓ OrdinaryDiffEqNonlinearSolve
3859.6 ms ✓ OrdinaryDiffEqStabilizedIRK
4120.6 ms ✓ OrdinaryDiffEqPDIRK
4835.3 ms ✓ OrdinaryDiffEqSDIRK
3541.2 ms ✓ OrdinaryDiffEqIMEXMultistep
3961.9 ms ✓ OrdinaryDiffEqFIRK
3302.9 ms ✓ OrdinaryDiffEqExponentialRK
8459.0 ms ✓ OrdinaryDiffEqBDF
45279.8 ms ✓ BoundaryValueDiffEqMIRK
16522.8 ms ✓ OrdinaryDiffEqDefault
5107.7 ms ✓ OrdinaryDiffEq
56868.6 ms ✓ BoundaryValueDiffEqFIRK
4447.0 ms ✓ DelayDiffEq
5782.5 ms ✓ BoundaryValueDiffEqShooting
6908.4 ms ✓ StochasticDiffEq
8007.4 ms ✓ BoundaryValueDiffEq
6748.5 ms ✓ DifferentialEquations
38 dependencies successfully precompiled in 123 seconds. 278 already precompiled.
Precompiling SciMLBaseMakieExt...
4202.7 ms ✓ SciMLBase → SciMLBaseMakieExt
1 dependency successfully precompiled in 5 seconds. 305 already precompiled.
Figure 2: Trajectory of the discontinuous system
We can also plot the trajectory in the state space, as in Fig. 3.
Figure 3: Trajectory of the discontinuous system in the state space
Now, how fast does the solution approach the origin?
Let’s use the 1-norm \|\bm x\|_1 = |x_1| + |x_2| to measure how far the trajectory is from the origin. We then ask:
\frac{\mathrm d}{\mathrm dt}\|\bm x\|_1 = ?
We avoid the troubles with nonsmoothness of the absolute value by consider each quadrant separately. Let’s start in the first (upper right) quadrant, that is, x_1>0 and x_2>0, and therefore |x_1| = x_1, \;|x_2| = x_2, and therefore
\frac{\mathrm d}{\mathrm dt}\|\bm x\|_1 = \dot x_1 + \dot x_2 = 1 - 3 = -2.
The situation is identical in the other quadrants. And, of course, undefined on the axes.
The conclusion is that the trajectory will hit the origin in finite time: for, say, x_1(0) = 1 and x_2(0) = 1 , the trajectory hits the origin at t=(|x_1(0)|+|x_2(0)|)/2 = 1. But with an infinite number of revolutions around the origin…
How will a standard algoritm for numerical simulation handle this? Let’s have a look at that.
Instead solving the above nonlinear equations with discontinuities, we introduce new variables u_1 and u_2 as the outputs of the \operatorname{sign} functions:
\begin{aligned}
{\color{blue} x_{1,k+1}} &= x_{1,k} + h (-{\color{blue}u_{1}} + 2 {\color{blue}u_{2}})\\
{\color{blue} x_{2,k+1}} &= x_{1,k} + h (-2{\color{blue}u_{1}} - {\color{blue}u_{2}}).
\end{aligned}
But now we have to enforce the relation between \bm u and \bm x_{k+1}. Recall the standard definition of the \operatorname{sign} function is
\operatorname{sign}(x) = \begin{cases}
1 & x>0\\
0 & x=0\\
-1 & x<0,
\end{cases}
but we change the definition to a set-valued function
\begin{cases}
\operatorname{sign}(x) = 1 & x>0\\
\operatorname{sign}(x) \in [-1,1] & x=0\\
\operatorname{sign}(x) = -1 & x<0.
\end{cases}
Accordingly, we set the relationship between \bm u and \bm x to
\begin{cases}
u_1 = 1 & x_1>0\\
u_1 \in [-1,1] & x_1=0\\
u_1 = -1 & x_1<0,
\end{cases}
and
\begin{cases}
u_2 = 1 & x_2>0\\
u_2 \in [-1,1] & x_2=0\\
u_2 = -1 & x_2<0.
\end{cases}
But these are mixed complementarity contraints we have defined previously!
\boxed{
\begin{aligned}
\begin{bmatrix}
{\color{blue} x_{1,k+1}}\\
{\color{blue} x_{1,k+1}}
\end{bmatrix}
&=
\begin{bmatrix}
x_{1,k}\\
x_{2,k}
\end{bmatrix} + h
\begin{bmatrix}
-1 & 2 \\
-2 & -1
\end{bmatrix}
\begin{bmatrix}
{\color{blue}u_{1}}\\
{\color{blue}u_{2}}
\end{bmatrix}\\
-1 \leq {\color{blue} u_1} \leq 1 \quad &\bot \quad -{\color{blue}x_{1,k+1}}\\
-1 \leq {\color{blue} u_2} \leq 1 \quad &\bot \quad -{\color{blue}x_{2,k+1}}.
\end{aligned}
}
9 possible combinations
There are now 9 possible combinations of the values of u_1 and u_2. Let’s explore some: x_{1,k+1} = x_{2,k+1} = 0, while u_1 \in [-1,1] and u_2 \in [-1,1]:
How does the set of states from which the next state is zero look like?
\begin{aligned}
-\begin{bmatrix}
-1 & 2 \\
-2 & -1
\end{bmatrix}^{-1}
\begin{bmatrix}
x_{1,k}\\
x_{2,k}
\end{bmatrix}
&= h
\begin{bmatrix}
{\color{blue}u_{1}}\\
{\color{blue}u_{2}}
\end{bmatrix}\\
-1 \leq {\color{blue} u_1} \leq 1, \quad -1 &\leq {\color{blue} u_2} \leq 1
\end{aligned}